## Packages

library(haven)
library(tidyr)
library(dplyr)
library(ggplot2)
library(ggpubr)
library(grid)
library(gridExtra)
library(xtable)
library(stargazer)
library(kableExtra)
library(RItools)
library(showtext)
library(texreg)
# library(screenreg)
library(estimatr)
library(tidylog)
library(pwr)
library(truncnorm)
library(InteractionPoweR)
library(ggplotify)
library(patchwork)
theme_set(theme_classic())

## Setting working directory 
# setwd("S:/Project/Delemappe/Samtykke")

#### Norwegian Citizen Panel #### 

# Importing dataset 
# citizens_full <- read_sav("Norsk medborgerpanel - runde 26 - v-104-L.sav")
# citizens_full <- citizens_full %>% select(responseid, r26_group, r26k2_lgcimx, 
                                          # r26k2_lgccax, r26k2_lgcimx_ran, r26P5_1, 
                                          # r26P1, r26P4_1, r26k2_bglrs, r26k2_bgurb)


# Filtering to respondents who received consent questions (r26_group = 4)
# citizens_full <- citizens_full %>% filter(r26_group == 4)

# Importing previous rounds for religion variable

# Importing Round 11
# citizens_r11 <- read_sav("Norsk medborgerpanel - runde 11 - v-102-L.sav")
# citizens_r11 <- citizens_r11 %>% select(responseid, r11bk41)
# citizens_r11 <- citizens_r11 %>% filter(r11bk41 < 9)

# Importing Round 14
# citizens_r14 <- read_sav("Norsk medborgerpanel - runde 14 - v-102-L.sav")
# citizens_r14 <- citizens_r14 %>% select(responseid, r14bk41)
# citizens_r14 <- citizens_r14 %>% filter(r14bk41 < 9)

# Importing Round 16
# citizens_r16 <- read_sav("Norsk medborgerpanel - runde 16 - v-102-L.sav")
# citizens_r16 <- citizens_r16 %>% select(responseid, r16bk41)
# citizens_r16 <- citizens_r16 %>% filter(r16bk41 < 9)

# Importing Round 18
# citizens_r18 <- read_sav("Norsk medborgerpanel - runde 18 - v-105-L.sav")
# citizens_r18 <- citizens_r18 %>% select(responseid, r18bk41)
# citizens_r18 <- citizens_r18 %>% filter(r18bk41 < 9)

# Importing Round 22
# citizens_r22 <- read_sav("Norsk medborgerpanel - runde 22 - v-102-L.sav")
# citizens_r22 <- citizens_r22 %>% select(responseid, r22bk41)
# citizens_r22 <- citizens_r22 %>% filter(r22bk41 < 9)

# Importing Round 25
# citizens_r25 <- read_sav("Norsk medborgerpanel - runde 25 - v-102-L.sav")
# citizens_r25 <- citizens_r25 %>% select(responseid, r25_bgrem)
# citizens_r25 <- citizens_r25 %>% filter(r25_bgrem < 9)

# Coding whether respondents are bureaucrats 

# Importing Round 23
# citizens_r23 <- read_sav("Norsk medborgerpanel - runde 23 - v-102-L.sav")
# citizens_r23 <- citizens_r23 %>% select(responseid, r23bk21_kodet)
# citizens_r23 <- citizens_r23 %>% 
  # select(responseid, r23bk21_kodet) %>% 
  # mutate(panel = ifelse(r23bk21_kodet == 1, "admin", "citizens"))

# Merging datasets
# citizens_full <- left_join(citizens_full, citizens_r11, by = "responseid")
# citizens_full <- left_join(citizens_full, citizens_r14, by = "responseid")
# citizens_full <- left_join(citizens_full, citizens_r16, by = "responseid")
# citizens_full <- left_join(citizens_full, citizens_r18, by = "responseid")
# citizens_full <- left_join(citizens_full, citizens_r22, by = "responseid")
# citizens_full <- left_join(citizens_full, citizens_r23, by = "responseid")
# citizens_full <- left_join(citizens_full, citizens_r25, by = "responseid")

# Removing responseid-variable for anonymity
# citizens_full <- citizens_full %>% select(!responseid)  

# Saving dataset
# citizens_full <- saveRDS(citizens_full, file = "citizens_full.rds")

# Loading dataset 
citizens_full <- readRDS("citizens_full.rds")

# Renaming and recoding variables 
citizens_full <- citizens_full %>% mutate(y1 = ifelse(r26k2_lgcimx == 97 | r26k2_lgcimx == 98, 
                                                       NA, (r26k2_lgcimx * -1) + 6), # Recoding direction of variable
                                           y1_std = (y1 - mean(y1, na.rm = TRUE)) / sd(y1, na.rm = TRUE), # Standardizing variable
                                           y2 = ifelse(r26k2_lgccax == 97 | r26k2_lgccax == 98, 
                                                       NA, (r26k2_lgccax * -1) + 8), # Recoding direction of variable
                                           y2_std = (y2 - mean(y2, na.rm = TRUE)) / sd(y2,  na.rm = TRUE), # Standardizing variable
                                           treatment = r26k2_lgcimx_ran, 
                                           age = ifelse(r26P5_1 == 6 | r26P5_1 == 7, 1, 0), # Young = 1, Old = 0
                                           age_con = r26P5_1, 
                                           age_cat = ifelse(r26P5_1 == 1, "1939 or earlier", 
                                                            ifelse(r26P5_1 == 2, "1940-1949", 
                                                                   ifelse(r26P5_1 == 3, "1950-1959",
                                                                          ifelse(r26P5_1 == 4, "1960-1969",
                                                                                 ifelse(r26P5_1 == 5, "1970-1979", 
                                                                                        ifelse(r26P5_1 == 6, "1980-1989", "1990 or later")))))),
                                           gender = ifelse(r26P1 == 1, 0, 1), # Man = 0, Woman = 1,
                                           education = ifelse(r26P4_1 > 2 & r26P4_1 < 97, 1, 
                                                              ifelse(r26P4_1 == 97, NA, 0)), # Recoding to binary variable
                                           lr = r26k2_bglrs, 
                                           urban = ifelse(r26k2_bgurb < 3, 1, 
                                                          ifelse(r26k2_bgurb > 3 & r26k2_bgurb < 97, 0, NA))) # Recoding to binary variable


# Coding religion variable
citizens_full <- citizens_full %>% mutate(religion = ifelse(r11bk41 < 9, r11bk41, NA),
                                          religion = ifelse(r14bk41 < 9 & is.na(religion), r14bk41, religion),
                                          religion = ifelse(r16bk41 < 9 & is.na(religion), r16bk41, religion),
                                          religion = ifelse(r18bk41 < 9 & is.na(religion), r18bk41, religion),
                                          religion = ifelse(r22bk41 < 9 & is.na(religion), r22bk41, religion),
                                          religion = ifelse(r25_bgrem < 9 & is.na(religion), r25_bgrem, religion), 
                                          religion = (religion*-1)+9) # Recoding direction of variable
table(citizens_full$religion, useNA = "always")
  

# Creating missing variables (response rate)
citizens_full$y1_RR <- ifelse(is.na(citizens_full$y1), 0, 1) 
citizens_full$y2_RR <- ifelse(is.na(citizens_full$y2), 0, 1)

table(citizens_full$treatment) # Checking that the number of respondents in each group is correct

# Pooling panels variable 
citizens_full <- citizens_full %>% mutate(panel = ifelse(is.na(panel), "citizens", panel))

### Subsetting data ###
# Creating baseline dataset (no information group)
citizens_baseline <- citizens_full %>% filter(treatment == 3)

# standardizing outcome variables with a mean of 0 and SD of 1
citizens_baseline <- citizens_baseline %>% mutate(y1_std = (y1 - mean(y1, na.rm = TRUE))/sd(y1, na.rm = TRUE), 
                                          y2_std = (y2 - mean(y2, na.rm = TRUE))/sd(y2, na.rm = TRUE))

# Creating main dataset (control and treatment group)
citizens <- citizens_full %>% filter(treatment == 1 | treatment == 2)

# standardizing outcome variables with a mean of 0 and SD of 1
citizens <- citizens %>% mutate(y1_std = (y1 - mean(y1, na.rm = TRUE))/sd(y1, na.rm = TRUE), 
                                                  y2_std = (y2 - mean(y2, na.rm = TRUE))/sd(y2, na.rm = TRUE))

# Recoding treatment variable to 0 = control and 1 = treatment 
citizens <- citizens %>% mutate(treatment = ifelse(treatment == 1, 0, 1))

# Creating dataset with respondents born after 1980 
citizens_young <- citizens %>% filter(age_cat == "1980-1989" | age_cat == "1990 or later")

# standardizing outcome variables with a mean of 0 and SD of 1
citizens_young <- citizens_young %>% mutate(y1_std = (y1 - mean(y1, na.rm = TRUE))/sd(y1, na.rm = TRUE), 
                                y2_std = (y2 - mean(y2, na.rm = TRUE))/sd(y2, na.rm = TRUE))

# Checking distributions between bureaucrats and normal citizens for pooled sample
table(citizens$panel, useNA = "always")

#### Norwegian Panel of Elected Representatives ####

# Importing dataset 
# elected_full <- read_sav("Representantpanelet - runde 9k2 - v-100-L.sav")
# elected_full <- elected_full %>% select(p9_group, p9k2_lgcimx, p9k2_lgccax, 
                                       # p9k2_lgcimx_ran, p9_P5_1, p9_P1, p9_P4_1,
                                       # p9k2_bglrs, p9k2_bgurb)
# elected_full <- saveRDS(elected_full, file = "elected_full.rds")

# Loading dataset
elected_full <- readRDS("elected_full.rds")

# Pooling panels variable 
elected_full$panel <- "elected"

# Renaming and recoding variables 
elected_full <- elected_full %>% mutate(y1 = ifelse(p9k2_lgcimx == 97 | p9k2_lgcimx == 98, 
                                                       NA, (p9k2_lgcimx * -1) + 6), # Recoding direction of variable
                                           y1_std = (y1 - mean(y1, na.rm = TRUE)) / sd(y1, na.rm = TRUE), # Standardizing variable
                                           y2 = ifelse(p9k2_lgccax == 97 | p9k2_lgccax == 98, 
                                                       NA, (p9k2_lgccax * -1) + 8), # Recoding direction of variable
                                           y2_std = (y2 - mean(y2, na.rm = TRUE)) / sd(y2,  na.rm = TRUE), # Standardizing variable
                                           treatment = p9k2_lgcimx_ran, 
                                           age = ifelse(p9_P5_1 == 5 | p9_P5_1 == 6, 1, 0), # Young = 1, Old = 0
                                           age_con = p9_P5_1, 
                                           age_cat =  ifelse(p9_P5_1 == 1, "1949 or earlier", 
                                                             ifelse(p9_P5_1 == 2, "1950-1959", 
                                                                    ifelse(p9_P5_1 == 3, "1960-1969",
                                                                           ifelse(p9_P5_1 == 4, "1970-1979",
                                                                                  ifelse(p9_P5_1 == 5, "1980-1989",
                                                                                         ifelse(p9_P5_1 == 6, "1990 or later", NA)))))),
                                           gender = ifelse(p9_P1 == 1, 0, 1), # Man = 0, Woman = 1
                                           education = ifelse(p9_P4_1 > 2 & p9_P4_1 < 97, 1, 
                                                              ifelse(p9_P4_1 == 97, NA, 0)), # Recoding to binary variable
                                           lr = p9k2_bglrs, 
                                           urban = ifelse(p9k2_bgurb < 3, 1,
                                                          ifelse(p9k2_bgurb > 3 & p9k2_bgurb < 97, 0, NA))) # Recoding to binary variable


# Filtering to respondents who received consent questions (r26_group = 4)
elected_full <- elected_full %>% filter(p9_group == 2)

# Creating missing variables (response rate)
elected_full$y1_RR <- ifelse(is.na(elected_full$y1), 0, 1) 
elected_full$y2_RR <- ifelse(is.na(elected_full$y2), 0, 1)

table(elected_full$treatment) # Checking that the number of respondents in each group is correct

### Subsetting data ###
# Creating baseline dataset (no information group)
elected_baseline <- elected_full %>% filter(treatment == 3)

# standardizing outcome variables with a mean of 0 and SD of 1
elected_baseline <- elected_baseline %>% mutate(y1_std = (y1 - mean(y1, na.rm = TRUE))/sd(y1, na.rm = TRUE), 
                              y2_std = (y2 - mean(y2, na.rm = TRUE))/sd(y2, na.rm = TRUE))


# Creating main dataset (control and treatment group)
elected <- elected_full %>% filter(treatment == 1 | treatment == 2)

# Recoding treatment variable to 0 = control and 1 = treatment 
elected <- elected %>% mutate(treatment = ifelse(treatment == 1, 0, 1))

# standardizing outcome variables with a mean of 0 and SD of 1
elected <- elected %>% mutate(y1_std = (y1 - mean(y1, na.rm = TRUE))/sd(y1, na.rm = TRUE), 
                                        y2_std = (y2 - mean(y2, na.rm = TRUE))/sd(y2, na.rm = TRUE))

#### Norwegian Panel of Public Administrators ####

# Importing dataset 
# administrators_full <- read_sav("Norsk forvaltningspanel - runde 3k2 - v-102-L.sav")
# administrators_full <- administrators_full %>% select(f3_group, f3k2_lgcimx, f3k2_lgccax, 
                                                   # f3k2_lgcimx_ran, f3_P5_1, f3_P1, f3_P4_1, 
                                                   # f3k2_bglrs, f3k2_bgurb)
# administrators_full <- saveRDS(administrators_full, file = "administrators_full.rds")

# Loading dataset
administrators_full <- readRDS("administrators_full.rds")

# Pooling panels variable 
administrators_full$panel <- "admin"

# Renaming and recoding variables 
administrators_full <- administrators_full %>% mutate(y1 = ifelse(f3k2_lgcimx == 97 | f3k2_lgcimx == 98, 
                                                    NA, (f3k2_lgcimx * -1) + 6), # Recoding direction of variable
                                        y1_std = (y1 - mean(y1, na.rm = TRUE)) / sd(y1, na.rm = TRUE), # Standardizing variable
                                        y2 = ifelse(f3k2_lgccax == 97 | f3k2_lgccax == 98, 
                                                    NA, (f3k2_lgccax * -1) + 8), # # Recoding direction of variable
                                        y2_std = (y2 - mean(y2, na.rm = TRUE)) / sd(y2,  na.rm = TRUE), # Standardizing variable
                                        treatment = f3k2_lgcimx_ran, 
                                        age = ifelse(f3_P5_1 == 4 | f3_P5_1 == 5, 1, 0), # Young = 1, Old = 0
                                        age_con = f3_P5_1, 
                                        age_cat =  ifelse(f3_P5_1 == 1, "1959 or earlier", 
                                                          ifelse(f3_P5_1 == 2, "1960-1969", 
                                                                 ifelse(f3_P5_1 == 3, "1970-1979",
                                                                        ifelse(f3_P5_1 == 4, "1980-1989",
                                                                               ifelse(f3_P5_1 == 5, "1990 or later", NA))))),
                                        gender = ifelse(f3_P1 == 1, 0, 1), # Man = 0, Woman = 1
                                        education = ifelse(f3_P4_1 > 2 & f3_P4_1 < 97, 1, 
                                                           ifelse(f3_P4_1 == 97, NA, 0)), # Recoding to binary variable
                                        lr = f3k2_bglrs, 
                                        urban = ifelse(f3k2_bgurb < 3, 1,
                                                       ifelse(f3k2_bgurb > 3 & f3k2_bgurb < 97, 0, NA))) # Recoding to binary variable


# Filtering to respondents who received consent questions (r26_group = 4)

administrators_full <- administrators_full %>% filter(f3_group == 2)

# Creating missing variables (response rate)
administrators_full$y1_RR <- ifelse(is.na(administrators_full$y1), 0, 1) 
administrators_full$y2_RR <- ifelse(is.na(administrators_full$y2), 0, 1)

table(administrators_full$treatment) # Checking that the number of respondents in each group is correct

# Creating panel variabel

### Subsetting data ###
# Creating baseline dataset (no information group)
administrators_baseline <- administrators_full %>% filter(treatment == 3)

# standardizing outcome variables with a mean of 0 and SD of 1
administrators_baseline <- administrators_baseline %>% mutate(y1_std = (y1 - mean(y1, na.rm = TRUE))/sd(y1, na.rm = TRUE), 
                                            y2_std = (y2 - mean(y2, na.rm = TRUE))/sd(y2, na.rm = TRUE))


# Creating main dataset (control and treatment group)
administrators <- administrators_full %>% filter(treatment == 1 | treatment == 2)

# Recoding treatment variable to 0 = control and 1 = treatment 
administrators <- administrators %>% mutate(treatment = ifelse(treatment == 1, 0, 1))

# standardizing outcome variables with a mean of 0 and SD of 1
administrators <- administrators %>% mutate(y1_std = (y1 - mean(y1, na.rm = TRUE))/sd(y1, na.rm = TRUE), 
                                                      y2_std = (y2 - mean(y2, na.rm = TRUE))/sd(y2, na.rm = TRUE))


## Figure 1: Distributions of dependent variables ---------------

### Distribution: Explicit consent ### 

showtext_opts(dpi = 300)

# Finding SD and mean for plot
sd_baseline_y1 <- sd(citizens_baseline$y1, na.rm = TRUE) %>% round(., digits = 2)
mean_baseline_y1 <- mean(citizens_baseline$y1, na.rm = TRUE) %>% round(., digits = 2)

# Plotting distribution
plot_baseline_Y1 <- citizens_baseline %>% 
  ggplot(.) + 
  geom_bar(aes(x = y1, y = after_stat(prop), group = 1), stat = "count", 
           fill = "lightgrey", color = "black") + 
  annotate("text", x = 0.84*mean_baseline_y1, y = 0.95*0.5, label = paste("Mean = ", mean_baseline_y1, 
                                                           sep = ""), size = 5) +
  annotate("text", x = 0.817*mean_baseline_y1, y = 0.87*0.5, label = paste("SD = ", sd_baseline_y1, 
                                                           sep = ""), size = 5) +
  geom_vline(xintercept = mean_baseline_y1, linetype = "dashed") + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5), 
                     breaks = seq(0, 0.5, by = 0.1),
                     labels = scales::percent_format()) +
  scale_x_continuous(breaks = seq(1, 5, by = 1), 
                     labels = c("1 \n Not at all \n important", "2 \n Not very \n important", 
                                "3 \n Somewhat \n important", "4 \n Important", 
                                "5 \n Very \n important")) + 
  labs(x = "Attitudes toward explicit consent", y = "", 
       title = "A. Explicit consent") + theme_classic() + 
  theme(axis.title.x = element_text(size = 20, margin = margin(t = 15, r = 0, b = 0, 
                                                               l = 0)),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 18),
        plot.title = element_text(vjust = 2, size = 20, face = "italic"),
        plot.title.position = "plot",
        plot.margin = margin(20, 20, 20, 20))

### Distribution: Consent scenario ###

# Finding SD and mean for plot
sd_baseline_y2 <- sd(citizens_baseline$y2, na.rm = TRUE) %>% round(., digits = 2)
mean_baseline_y2 <- mean(citizens_baseline$y2, na.rm = TRUE) %>% round(., digits = 2)

# Plotting baseline distribution
plot_baseline_Y2 <- citizens_baseline %>%
  ggplot(.) + 
  geom_bar(aes(x = y2, y = after_stat(prop), group = 1), stat = "count", 
           fill = "lightgrey", color = "black") + 
  annotate("text", x = 0.805*mean_baseline_y2, y = 0.95*0.4, label = paste("Mean = ", mean_baseline_y2, 
                                                           sep = ""), size = 5) +
  annotate("text", x = 0.778*mean_baseline_y2, y = 0.87*0.4, label = paste("SD = ", sd_baseline_y2, 
                                                           sep = ""), size = 5) +
  geom_vline(xintercept = mean_baseline_y2, linetype = "dashed") + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, 0.4),
                     breaks = seq(0, 0.4, by = 0.1),
                     labels = scales::percent_format()) +
  scale_x_continuous(breaks = seq(1, 7, by = 1), 
                     labels = c("1 \n Completely \n disagree", "2 \n Disagree", 
                                "3 \n Somewhat \n disagree",
                                "4 \n Neither \n / or", 
                                "5 \n Somewhat \n agree", 
                                "6 \n Agree", "7 \n Completely \n agree")) + 
  labs(x = "Attitudes toward consent scenario", y = "", 
       title = "B. Consent scenario") + theme_classic() + 
  theme(axis.title.x = element_text(size = 20, margin = margin(t = 15, r = 0, b = 0, 
                                                                l = 0)),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 18),
        plot.title = element_text(vjust = 2, size = 20, face = "italic"),
        plot.title.position = "plot",
        plot.margin = margin(20, 20, 20, 20))

### Combining plots ###

grid.arrange(plot_baseline_Y1, plot_baseline_Y2, nrow = 1)
distribution_baselines <- arrangeGrob(plot_baseline_Y1, plot_baseline_Y2, nrow = 1)
ggsave(height = 5, width = 15, filename = "Fig_1_distribution_baselines.png", 
       distribution_baselines)

### Descriptives ###

prop.table(table(citizens_baseline$y1))
prop.table(table(citizens_baseline$y2))


## Tables L.1 and L.2: Baseline attitudes correlated with pre-treatment covariates ----

### Regression: Baseline (no information group): 1939 or earlier as reference ###

# Model with covariates: Explicit consent 
baseline_y1 <- lm_robust(y1_std ~ gender + age_cat + education + urban + 
                         religion +  lr, data = citizens_baseline)
summary(baseline_y1)

# Model with covariates: Consent scenario 
baseline_y2 <- lm_robust(y2_std ~ gender + age_cat + education + urban + 
                           religion + lr, data = citizens_baseline)
summary(baseline_y2)

# Table for paper
screenreg(list(baseline_y1, baseline_y2), 
          include.ci = FALSE, 
          custom.coef.map = list("(Intercept)" = "Constant", 
                                 "gender" = "Woman", 
                                 "age_cat1939 or earlier" = "Born 1949 or earlier", 
                                 "age_cat1940-1949" = "Born 1940-1949", 
                                 "age_cat1950-1959" = "Born 1950-1959",
                                 "age_cat1960-1969" = "Born 1960-1969",
                                 "age_cat1970-1979" = "Born 1970-1979", 
                                 "age_cat1980-1989" = "Born 1980-1989",
                                 "age_cat1990 or later" = "Born 1990 or later",
                                 "education" = "Education", "urban" = "Urban", 
                                 "religion" = "Religiosity (higher values = more religious)", 
                                 "lr" = "Left-right orientation (higher values = right)"))

texreg(list(baseline_y1, baseline_y2), 
       file = "Table_L1_baseline_attitudes.tex",
       include.ci = FALSE, 
       include.rmse = FALSE, 
       caption = "Baseline attitudes correlated with pre-treatment covariates",
       label = "tab:baseline_correlates",
       caption.above = TRUE, 
       custom.header = list("\\emph{Dependent variable}" = 1:2),
       custom.model.name = c("Explicit consent", "Consent scenario"),
       custom.coef.map = list("(Intercept)" = "Constant", 
                              "gender" = "Woman", 
                              "age_cat1939 or earlier" = "Born 1949 or earlier", 
                              "age_cat1940-1949" = "Born 1940-1949", 
                              "age_cat1950-1959" = "Born 1950-1959",
                              "age_cat1960-1969" = "Born 1960-1969",
                              "age_cat1970-1979" = "Born 1970-1979", 
                              "age_cat1980-1989" = "Born 1980-1989",
                              "age_cat1990 or later" = "Born 1990 or later",
                              "education" = "Higher education", "urban" = "Urban residence", 
                              "religion" = "Religiosity (higher values = more religious)", 
                              "lr" = "Left-right orientation (higher values = right)"), 
       threeparttable = TRUE,
       custom.note = "\n\\item \\footnotesize \\textit{Note:} OLS regression models with standard errors in parentheses. Sample consists of respondents 
       \\item \\footnotesize from the Norwegian Citizen Panel. Man is the reference category for gender. Born in 1939 or
       \\item \\footnotesize earlier is the reference category for age. Less than higher education is the reference category for 
       \\item \\footnotesize education. Rural residence is the reference category for residence. DVs are standardized
       \\item \\footnotesize with a mean of 0 and a SD of 1. Significant at %stars.")

### Regression: Baseline (no information group): 1990 or later as reference ###

citizens_baseline <- citizens_baseline %>% mutate(age_cat_1990 = factor(age_cat, levels = c("1990 or later", 
                                                                                            "1980-1989", "1970-1979",
                                                                                            "1960-1969", "1950-1959", 
                                                                                            "1940-1949", "1939 or earlier")))

# Model with covariates: Explicit consent 
baseline_y1_1990 <- lm_robust(y1_std ~ gender + age_cat_1990 + education + urban + 
                           religion +  lr, data = citizens_baseline)
summary(baseline_y1_1990)

# Model with covariates: Consent scenario 
baseline_y2_1990 <- lm_robust(y2_std ~ gender + age_cat_1990 + education + urban + 
                           religion + lr, data = citizens_baseline)
summary(baseline_y2_1990)

# Table for paper
screenreg(list(baseline_y1_1990, baseline_y2_1990), 
          include.ci = FALSE, 
          custom.coef.map = list("(Intercept)" = "Constant", 
                                 "gender" = "Woman", 
                                 "age_cat_19901939 or earlier" = "Born 1939 or earlier", 
                                 "age_cat_19901940-1949" = "Born 1940-1949", 
                                 "age_cat_19901950-1959" = "Born 1950-1959",
                                 "age_cat_19901960-1969" = "Born 1960-1969",
                                 "age_cat_19901970-1979" = "Born 1970-1979", 
                                 "age_cat_19901980-1989" = "Born 1980-1989",
                                 "age_cat_19901990 or later" = "Born 1990 or later",
                                 "education" = "Higher education", "urban" = "Urban residence", 
                                 "religion" = "Religiosity (higher values = more religious)", 
                                 "lr" = "Left-right orientation (higher values = right)"))

texreg(list(baseline_y1_1990, baseline_y2_1990), 
       file = "Table_L2_baseline_attitudes_1990.tex",
       include.ci = FALSE, 
       include.rmse = FALSE, 
       caption = "Baseline attitudes correlated with pre-treatment covariates",
       label = "tab:baseline_correlates",
       caption.above = TRUE, 
       custom.header = list("\\emph{Dependent variable}" = 1:2),
       custom.model.name = c("Explicit consent", "Consent scenario"),
       custom.coef.map = list("(Intercept)" = "Constant", 
                              "gender" = "Woman", 
                              "age_cat_19901939 or earlier" = "Born 1939 or earlier", 
                              "age_cat_19901940-1949" = "Born 1940-1949", 
                              "age_cat_19901950-1959" = "Born 1950-1959",
                              "age_cat_19901960-1969" = "Born 1960-1969",
                              "age_cat_19901970-1979" = "Born 1970-1979", 
                              "age_cat_19901980-1989" = "Born 1980-1989",
                              "age_cat_19901990 or later" = "Born 1990 or later",
                              "education" = "Higher education", "urban" = "Urban residence", 
                              "religion" = "Religiosity (higher values = more religious)", 
                              "lr" = "Left-right orientation (higher values = right)"), 
       threeparttable = TRUE,
       custom.note = "\n\\item \\footnotesize \\textit{Note:} OLS regression models with standard errors in parentheses. Sample consists of respondents 
       \\item \\footnotesize from the Norwegian Citizen Panel. Man is the reference category for gender. Born in 1939 or
       \\item \\footnotesize earlier is the reference category for age. Less than higher education is the reference category for 
       \\item \\footnotesize education. Rural residence is the reference category for residence. DVs are standardized
       \\item \\footnotesize with a mean of 0 and a SD of 1. Significant at %stars.")


## Tables J.1-3: Balance tests ---------------------------------------------

### Balance test: Norwegian Citizen Panel ###

# Control group and treatment group
bt_citizens <- balanceTest(treatment ~ gender + age + education +
                             urban + religion + lr, data = citizens)

bt_citizens_tab <- xtable(bt_citizens,
                          caption = "Balance test for the Citizen Panel: control group versus treatment group")

output_bt_citizens <- print(bt_citizens_tab, 
                            include.rownames = TRUE,
                            hline.after = c(-1, 0, nrow(data)),
                            size = "small", 
                            print.results = FALSE)

chi_square_citizens <-  paste("\\\\hline\n\\\\multicolumn{7}{l}{\\\\emph{Overall test}} \\\\\\\\\n\\\\hline\nChi-square &", 
                              round(bt_citizens$overall$chisquare, digits = 2), 
                              "& & & & & \\\\\\\\\nDf &", bt_citizens$overall$df ,
                              "& & & & & \\\\\\\\\nP-value &", 
                              round(bt_citizens$overall$p.value, digits = 2), 
                              "& & & & & \\\\\\\\ \\\\hline\n\\\\end{tabular}", 
                              sep = "")

output_bt_citizens <- gsub("\\\\end\\{tabular\\}", chi_square_citizens, 
                           output_bt_citizens)

cat(output_bt_citizens, 
    file = "Table_J1_balance_citizen_panel.tex")

### Balance test: Panel of Public Administrators ###

# Control group and treatment group
bt_administrators <- balanceTest(treatment ~ gender + age + education + 
                                   urban + lr, data = administrators)

bt_administrators_tab <- xtable(bt_administrators, 
                                caption = "Balance test for the Panel of Public Administrators: control group versus treatment group")

output_bt_administrators <- print(bt_administrators_tab, 
                                  include.rownames = TRUE,
                                  hline.after = c(-1, 0, nrow(data)),
                                  size = "small", 
                                  print.results = FALSE)

chi_square_administrators <-  paste("\\\\hline\n\\\\multicolumn{7}{l}{\\\\emph{Overall test}} \\\\\\\\\n\\\\hline\nChi-square &", 
                                    round(bt_administrators$overall$chisquare, digits = 2), 
                                    "& & & & & \\\\\\\\\nDf &", bt_administrators$overall$df ,
                                    "& & & & & \\\\\\\\\nP-value &", 
                                    round(bt_administrators$overall$p.value, digits = 2), 
                                    "& & & & & \\\\\\\\ \\\\hline\n\\\\end{tabular}", 
                                    sep = "")
output_bt_administrators <- gsub("\\\\end\\{tabular\\}", chi_square_administrators, 
                                 output_bt_administrators)

cat(output_bt_administrators, file = "Table_J2_balance_pub_admin.tex")

### Balance test: Panel of Elected Representatives ###

# Control group and treatment group
bt_elected <- balanceTest(treatment ~ gender + age + education + 
                            urban + lr, data = elected)

bt_elected_tab <- xtable(bt_elected, 
                         caption = "Balance test for the Panel of Elected Representatives: control group versus treatment group")

output_bt_elected <- print(bt_elected_tab, 
                           include.rownames = TRUE,
                           hline.after = c(-1, 0, nrow(data)),
                           size = "small", 
                           print.results = FALSE)

chi_square_elected <-  paste("\\\\hline\n\\\\multicolumn{7}{l}{\\\\emph{Overall test}} \\\\\\\\\n\\\\hline\nChi-square &", 
                             round(bt_elected$overall$chisquare, digits = 2), 
                             "& & & & & \\\\\\\\\nDf &", bt_elected$overall$df ,
                             "& & & & & \\\\\\\\\nP-value &", 
                             round(bt_elected$overall$p.value, digits = 2), 
                             "& & & & & \\\\\\\\ \\\\hline\n\\\\end{tabular}", 
                             sep = "")
output_bt_elected <- gsub("\\\\end\\{tabular\\}", chi_square_elected, 
                          output_bt_elected)

cat(output_bt_elected, file = "Table_J3_balance_elect_repr.tex")


table(citizens$education)

## Tables K.1-2: Differential attrition ------------------------------------

### Table K.1: Attrition rate: no information group vs. control group ###

# Citizen panel

citizens_BC <- citizens_full %>% filter(treatment == 1 | treatment == 3) %>% 
  mutate(treatment = factor(treatment)) 

citizens_BC_y1 <- t.test(y1_RR ~ treatment, data = citizens_BC) # Explicit consent 
citizens_BC_y2 <- t.test(y2_RR ~ treatment, data = citizens_BC) # Consent scenario

# Panel of Public Administrators:

administrators_BC <- administrators_full %>% filter(treatment == 1 | treatment == 3) %>% 
  mutate(treatment = factor(treatment)) 

administrators_BC_y1 <- t.test(y1_RR ~ treatment, data = administrators_BC) # Explicit consent 
administrators_BC_y2 <- t.test(y2_RR ~ treatment, data = administrators_BC) # Consent scenario

# Panel of Elected Officials

elected_BC <- elected_full %>% filter(treatment == 1 | treatment == 3) %>% 
  mutate(treatment = factor(treatment)) 

elected_BC_y1 <- t.test(y1_RR ~ treatment, data = elected_BC) # Explicit consent 
elected_BC_y2 <- t.test(y2_RR ~ treatment, data = elected_BC) # Consent scenario

# Creating table 

attrition_labels <- c("citizens_y1", "citizens_y2", "administrators_y1", "administrators_y2",
                      "elected_y1", "elected_y2")
y_labels <- c("Explicit consent", "Consent scenario", 
              "Explicit consent", "Consent scenario", 
              "Explicit consent", "Consent scenario")

baseline_attrition <- c((1 - citizens_BC_y1$estimate[2])*100, 
                        (1 - citizens_BC_y2$estimate[2])*100, 
                        (1 - administrators_BC_y1$estimate[2])*100,
                        (1 - administrators_BC_y2$estimate[2])*100, 
                        (1 - elected_BC_y1$estimate[2])*100, 
                        (1 - elected_BC_y2$estimate[2])*100)
baseline_attrition <- baseline_attrition %>% unname() %>% round(., digits = 2)

control_attrition <- c((1 - citizens_BC_y1$estimate[1])*100, 
                       (1 - citizens_BC_y2$estimate[1])*100, 
                       (1 - administrators_BC_y1$estimate[1])*100,
                       (1 - administrators_BC_y2$estimate[1])*100, 
                       (1 - elected_BC_y1$estimate[1])*100, 
                       (1 - elected_BC_y2$estimate[1])*100)
control_attrition <- control_attrition %>% unname() %>% round(., digits = 2)

BC_attrition <- data.frame(baseline_attrition, control_attrition) # Added to get correct difference in attrition

BC_difference <- BC_attrition$baseline_attrition - BC_attrition$control_attrition  # Modified to get correct difference in attrition

BC_p.value <- c(citizens_BC_y1$p.value, citizens_BC_y2$p.value, 
                administrators_BC_y1$p.value, administrators_BC_y2$p.value, 
                elected_BC_y1$p.value,elected_BC_y2$p.value) %>% round(., digits = 2)

BC_data <- data.frame(attrition_labels, y_labels, baseline_attrition, control_attrition,
                      BC_difference, BC_p.value)

BC_data[1:6, 2:6] %>% kbl(format="latex",
                          col.names = c("","B","C","B - C","p-value"),
                          align="r", 
                          booktabs = T) %>% kable_styling() %>%
  pack_rows("Norwegian Citizen Panel", 1, 2, bold = F, italic = T) %>%
  pack_rows("Panel of Public Administrators", 3, 4, bold = F, italic = T, hline_before = TRUE) %>%
  pack_rows("Panel of Elected Representatives", 5, 6, bold = F, italic = T, hline_before = TRUE) %>%
  save_kable("Table_K1_attrition_baseline_control.tex")


### Table K.2: Attrition rate: no information group vs. treatment group ###

# Citizen Panel:
citizens_BT <- citizens_full %>% filter(treatment == 2 | treatment == 3)

citizens_BT_y1 <- t.test(y1_RR ~ treatment, data = citizens_BT) # Explicit consent 
citizens_BT_y2 <- t.test(y2_RR ~ treatment, data = citizens_BT) # Consent scenario

# Panel of Public Administrators
administrators_BT <- administrators_full %>% filter(treatment == 2 | treatment == 3)

administrators_BT_y1 <- t.test(y1_RR ~ treatment, data = administrators_BT) # Explicit consent 
administrators_BT_y2 <- t.test(y2_RR ~ treatment, data = administrators_BT) # Consent scenario

# Panel of Elected Officials
elected_BT <- elected_full %>% filter(treatment == 2 | treatment == 3)

# elected_BT_y1 <- t.test(y1_RR ~ treatment, data = elected_BT) # Explicit consent 
# No missing on y1 for elected representatives
table(elected_BT$y2_RR, elected_BT$treatment, useNA = "always")

elected_BT_y2 <- t.test(y2_RR ~ treatment, data = elected_BT) # Consent scenario

# Creating table
treatment_attrition <- c((1 - citizens_BT_y1$estimate[1])*100, 
                         (1 - citizens_BT_y2$estimate[1])*100, 
                         (1 - administrators_BT_y1$estimate[1])*100,
                         (1 - administrators_BT_y2$estimate[1])*100, 
                          0, # No missingness for y1 and elected
                         (1 - elected_BT_y2$estimate[1])*100)
treatment_attrition <- treatment_attrition %>% unname() %>% round(., digits = 2)

BT_attrition <- data.frame(baseline_attrition, treatment_attrition) # Added to get correct difference in attrition

BT_difference <- BT_attrition$baseline_attrition - BT_attrition$treatment_attrition # Modified to get correct difference in attrition

BT_p.value <- c(citizens_BT_y1$p.value, citizens_BT_y2$p.value, 
                administrators_BT_y1$p.value, administrators_BT_y2$p.value, 
                NA, # No missingness for y1 and elected
                elected_BT_y2$p.value) %>% round(., digits = 2)

BT_data <- data.frame(attrition_labels, y_labels, baseline_attrition, treatment_attrition,
                      BT_difference, BT_p.value)

BT_data[1:6, 2:6] %>% kbl(format="latex",
                          col.names = c("","B","T","B - T","p-value"),
                          align="r", 
                          booktabs = T) %>% kable_styling() %>%
  pack_rows("Norwegian Citizen Panel", 1, 2, bold = F, italic = T) %>%
  pack_rows("Panel of Public Administrators", 3, 4, bold = F, italic = T, hline_before = TRUE) %>%
  pack_rows("Panel of Elected Representatives", 5, 6, bold = F, italic = T, hline_before = TRUE) %>%
  save_kable("Table_K2_attrition_baseline_treatment.tex")


## Figure 2 and Table M.1: Average treatment effect --------

### Regression models ###

# Explicit consent as dependent variable
h1_y1 <- lm_robust(y1_std ~ treatment, data = citizens)
summary(h1_y1)

tapply(citizens$y1_std, citizens$treatment, mean, na.rm = TRUE)

# Consent scenario as dependent variable
h1_y2 <- lm_robust(y2_std ~ treatment, data = citizens)
summary(h1_y2)

### Table M.1: Testing H_consent: average treatment effect ###

screenreg(list(h1_y1, h1_y2), include.ci = FALSE, 
          custom.coef.map = list("(Intercept)" = "Constant", "treatment" = "Treatment"),
          custom.model.names = c("Explicit consent", "Consent scenario"))

texreg(list(h1_y1, h1_y2), 
       include.ci = FALSE, 
       include.rmse = FALSE, 
       file = "Table_M1_average_treatment_effect.tex",
       caption = "Testing $H_{consent}$: average treatment effect",
       label = "tab:ate_overall",
       caption.above = TRUE, 
       custom.header = list("Dependent variable" = 1:2),
       custom.model.name = c("Explicit consent", "Consent scenario"),
       custom.coef.map = list("(Intercept)" = "Constant", "treatment" = "Treatment"), 
       threeparttable = TRUE,
       custom.note = "\n\\item \\footnotesize \\textit{Note:} OLS regression models with standard errors in parentheses. Sample consists of
       \\item \\footnotesize respondents from the Norwegian Citizen Panel. DVs are standardized with a mean of
       \\item \\footnotesize  0 and a SD of 1. Significant at %stars."
)

### Coefficient plot ###

# Explicit consent as dependent variable
h1_y1_lh <- lh_robust(y1_std ~ treatment,
                      data = citizens,
                      linear_hypothesis = "treatment = 0")
tidy(h1_y1_lh)

dt.p <- bind_rows(tidy(h1_y1_lh) %>% filter(term == "treatment = 0"))

v.pval <- tidy(h1_y1_lh) %>% 
  filter(term == "treatment") %>% 
  dplyr::select(p.value) %>% 
  pull()

dt.p <- dt.p %>% 
  mutate(diff_p = v.pval,
         xmin = 1,
         xmax = 2, 
         ypos = max(conf.high), 
         p.value_bh = p.value * 2)

plot_h1_y1 <- dt.p %>% 
  mutate(term = ifelse(term == "treatment = 0", "Overall \n sample", NA)) %>%
  ggplot(aes(y = estimate, x = term, ymin = conf.low,
             ymax = conf.high)) + 
  geom_point(size = 4) +
  geom_errorbar(width = 0.1) + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(y = "Treatment effect (with 95% CIs)",
       x = "",
       title = "A. Explicit consent") + 
  scale_y_continuous(expand= c(0,0), limits = c(-0.5, 0.1),
                     breaks = seq(-0.5, 0.1, 0.1)) +
  coord_flip() +
  theme_minimal() +  
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0, size = 24, face = "italic"),
        plot.subtitle = element_text(hjust = 0, size = 15, face = "italic"),
        plot.title.position = "plot",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(),
        axis.line.x = element_line(linewidth = .3),
        axis.title = element_text(size = 22),
        axis.text = element_text(size = 20),
        axis.text.x = element_text(vjust = -2),
        axis.title.x = element_text(vjust = -3, margin = margin(t = 10, r = 0, b = 0, 
                                                              l = 0)),
        plot.margin = margin(16, 20, 18, 18), 
        aspect.ratio = 0.3) 

# Consent scenario as dependent variable
h1_y2_lh <- lh_robust(y2_std ~ treatment,
                      data = citizens,
                      linear_hypothesis = "treatment = 0")
tidy(h1_y2_lh)

dt.p <- bind_rows(tidy(h1_y2_lh) %>% filter(term == "treatment = 0"))

v.pval <- tidy(h1_y2_lh) %>% 
  filter(term == "treatment") %>% 
  dplyr::select(p.value) %>% 
  pull()

dt.p <- dt.p %>% 
  mutate(diff_p = v.pval,
         xmin = 1,
         xmax = 2, 
         ypos = max(conf.high), 
         p.value_bh = p.value * 2)

plot_h1_y2 <- dt.p %>% 
  mutate(term = ifelse(term == "treatment = 0", "Overall \n sample", NA)) %>%
  ggplot(aes(y = estimate, x = term, ymin = conf.low,
             ymax = conf.high)) + 
  geom_point(size = 4) +
  geom_errorbar(width = 0.1) + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(y = "Treatment effect (with 95% CIs)",
       x = "",
       title = "B. Consent scenario") + 
  scale_y_continuous(expand= c(0,0), limits = c(-0.5, 0.1),
                     breaks = seq(-0.5, 0.1, 0.1)) +
  coord_flip() +
  theme_minimal() +  
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0, size = 24, face = "italic"),
        plot.subtitle = element_text(hjust = 0, size = 15, face = "italic"),
        plot.title.position = "plot",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(),
        axis.line.x = element_line(linewidth = .3),
        axis.title = element_text(size = 22),
        axis.text = element_text(size = 20),
        axis.text.x = element_text(vjust = -2),
        axis.title.x = element_text(vjust = -3, margin = margin(t = 10, r = 0, b = 0, 
                                                                l = 0)),
        plot.margin = margin(16, 20, 18, 18), 
        aspect.ratio = 0.3)

### Combining plots ###

grid.arrange(plot_h1_y1, plot_h1_y2, nrow = 1)
ate_overall <- arrangeGrob(plot_h1_y1, plot_h1_y2, nrow = 1)
ggsave(height = 5, width = 15, filename = "Fig_2_avg_treatment_eff.png", 
       ate_overall)

## Tables O.1-2: Robustness tests: analysis w/reduced sample -----------------

### Table O.1: Removing extreme values ###

# Removing the extreme values
citizens_extreme_y1 <- citizens %>% filter(y1 > 1.5) 

# Rerunning main regression models
# Explicit consent as dependent variable
h1_y1_extreme <- lm_robust(y1_std ~ treatment, data = citizens_extreme_y1)
summary(h1_y1_extreme)

# Consent scenario as dependent variable
h1_y2_extreme <- lm_robust(y2_std ~ treatment, data = citizens_extreme_y1)
summary(h1_y2_extreme)

# Table output
screenreg(list(h1_y1_extreme, h1_y2_extreme), include.ci = FALSE, 
          custom.coef.map = list("(Intercept)" = "Constant", "treatment" = "Treatment"),
          custom.model.names = c("Explicit consent", "Consent scenario"))

texreg(list(h1_y1_extreme, h1_y2_extreme), 
       include.ci = FALSE, 
       include.rmse = FALSE, 
       file = "Table_O1_removing_extreme_values.tex",
       caption = "Testing $H_{consent}$: removing extreme values",
       label = "tab:ate_extreme",
       caption.above = TRUE, 
       custom.header = list("Dependent variable" = 1:2),
       custom.coef.map = list("(Intercept)" = "Constant", "treatment" = "Treatment"),
       custom.model.names = c("Explicit consent", "Consent scenario"), 
       threeparttable = TRUE,
       custom.note = "\n\\item \\footnotesize \\textit{Note:} OLS regression models with standard errors in parentheses. Sample consists of
       \\item \\footnotesize respondents from the Norwegian Citizen Panel excluding those answering ``Not at
       \\item \\footnotesize all important on the explicit consent variable. DVs are standardized with a mean of
       \\item \\footnotesize 0 and a SD of 1. Significant at %stars.")

### Table O.2: removing individuals with large discrepancies

# Excluding those with at least a 4-point discrepancy between
# their answers to the two dependent variables (explicit
# consent and consent scenario)
citizens <- citizens %>% mutate(difference = y1 - y2)

citizens_difference <- citizens %>% filter(difference <= 4)

# Rerunning main regression models
# Explicit consent as dependent variable 
h1_y1_difference <- lm_robust(y1_std ~ treatment, data = citizens_difference)
summary(h1_y1_difference)

# Consent scenario as dependent variable
h1_y2_difference <- lm_robust(y2_std ~ treatment, data = citizens_difference)
summary(h1_y2_difference)

screenreg(list(h1_y1_difference, h1_y2_difference), include.ci = FALSE, 
          custom.coef.map = list("(Intercept)" = "Constant", "treatment" = "Treatment"),
          custom.model.names = c("Explicit consent", "Consent scenario"))

texreg(list(h1_y1_difference, h1_y2_difference), 
       include.ci = FALSE, 
       include.rmse = FALSE, 
       file = "Table_O2_removing_large_discrepancies.tex",
       caption = "Testing $H_{consent}$: removing individuals with large discrepancies",
       label = "tab:ate_difference",
       caption.above = TRUE, 
       custom.header = list("Dependent variable" = 1:2),
       custom.coef.map = list("(Intercept)" = "Constant", "treatment" = "Treatment"),
       custom.model.names = c("Explicit consent", "Consent scenario"),
       threeparttable = TRUE,
       custom.note = "\n\\item \\footnotesize \\textit{Note:} OLS regression models with standard errors in parentheses. Sample consists of
       \\item \\footnotesize respondents from the Norwegian Citizen Panel excluding those with at least a 4-point
       \\item \\footnotesize discrepancy between their answers to our two questions. DVs are standardized with a 
       \\item \\footnotesize mean of 0 and a SD of 1. Significant at %stars.")


## Figure 3 and Tables M.2-3: Exploring heterog. treatment eff. --------

### Table M.2: Heterogeneous effects for the young and old ###

# Explicit consent
h2a_y1 <- lm_robust(y1_std ~ treatment*age, data = citizens)
summary(h2a_y1)

# Consent scenario
h2a_y2 <- lm_robust(y2_std ~ treatment*age, data = citizens)
summary(h2a_y2)

# Table output
screenreg(list(h2a_y1, h2a_y2), include.ci = FALSE, 
          custom.coef.map = list("(Intercept)" = "Constant", "treatment" = "Treatment", 
                                 "age" = "Born in/after 1980", 
                                 "treatment:age" = "Treatment*Born in/after 1980"))

texreg(list(h2a_y1, h2a_y2), 
       include.ci = FALSE, 
       include.rmse = FALSE, 
       file = "Table_M2_heterogenous_age.tex",
       caption = "Exploring heterogeneous treatment effects for the young (born in/after 1980) and the old (born before 1980)",
       label = "tab:cate_age",
       caption.above = TRUE, 
       custom.header = list("Dependent variable" = 1:2),
       custom.model.name = c("Explicit consent", "Consent scenario"),
       custom.coef.map = list("(Intercept)" = "Constant", "treatment" = "Treatment", 
                              "age" = "Born in/after 1980", 
                              "treatment:age" = "Treatment*Born in/after 1980"), 
       threeparttable = TRUE,
       custom.note = "\n\\item \\footnotesize \\textit{Note:} OLS regression models with standard errors in parentheses. Sample consists of
       \\item \\footnotesize respondents from the Norwegian Citizen Panel. DVs are standardized with a mean of
       \\item \\footnotesize  0 and a SD of 1. Significant at %stars.")

### Table M.3: Heterogeneous effects for young women and men ###

# Explicit consent as dependent variable
h2b_y1 <- lm_robust(y1_std ~ treatment*gender, data = citizens_young)
summary(h2b_y1)

# Consent scenario as dependent variable
h2b_y2 <- lm_robust(y2_std ~ treatment*gender, data = citizens_young)
summary(h2b_y2)

# Table for paper
screenreg(list(h2b_y1, h2b_y2), include.ci = FALSE, 
          custom.coef.map = list("(Intercept)" = "Constant", "treatment" = "Treatment", 
                                 "gender" = "Woman",
                                 "treatment:gender" = "Treatment*Woman"))

texreg(list(h2b_y1, h2b_y2), 
       include.ci = FALSE, 
       include.rmse = FALSE, 
       file = "Table_M3_heterogenous_gender.tex",
       caption = "Exploring heterogeneous treatment effect for young women and young men",
       label = "tab:cate_gender",
       caption.above = TRUE, 
       custom.header = list("Dependent variable" = 1:2),
       custom.model.name = c("Explicit consent", "Consent scenario"),
       custom.coef.map = list("(Intercept)" = "Constant", "treatment" = "Treatment", 
                              "gender" = "Woman",
                              "treatment:gender" = "Treatment*Woman"), 
       threeparttable = TRUE,
       custom.note = "\n\\item \\footnotesize \\textit{Note:} OLS regression models with standard errors in parentheses. Sample consists of
       \\item \\footnotesize respondents from the Norwegian Citizen Panel born in/after 1980. DVs are standardized
       \\item \\footnotesize with a mean of 0 and a SD of 1. Significant at %stars.")


### Figure 3: Conditional average treatment effect by age and gender ###

# Coefficient plots for age heterogeneous effects (Table M.2)

# Explicit consent
h2a_y1_lh <- lh_robust(y1_std ~ treatment*age,
                       data = citizens,
                       linear_hypothesis = "treatment + treatment:age = 0")
tidy(h2a_y1_lh)

dt.p <- bind_rows(tidy(h2a_y1_lh) %>% filter(term=="treatment"),
                  tidy(h2a_y1_lh) %>% filter(term=="treatment + treatment:age = 0"))

v.pval <- tidy(h2a_y1_lh) %>% 
  filter(term == "treatment:age") %>% 
  dplyr::select(p.value) %>% 
  pull()

dt.p <- dt.p %>% 
  mutate(diff_p = v.pval,
         xmin = 1,
         xmax = 2, 
         ypos = max(conf.high))

plot_h2a_y1 <- dt.p %>% 
  mutate(label = if_else(term == "treatment", "Old", "Young")) %>% 
  ggplot(aes(y = estimate, x = label, ymin = conf.low,
             ymax = conf.high)) + 
  geom_point(size = 4) +
  geom_errorbar(width = 0.1) + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_bracket(data=dt.p %>% slice(1),
               aes(xmin = xmin,
                   xmax = xmax,
                   y.position = ypos,
                   label=if_else(diff_p<.001,
                                 "p < 0.001",
                                 paste0("p = ",(round(diff_p,2))))),
               coord.flip = TRUE,
               label.size = 6,
               color="#566573",
               tip.length = .02,
               vjust = -.7,
               bracket.nudge.y = 0.05,
               inherit.aes = FALSE) +
  labs(y = "Treatment effect (with 95% CIs)",
       x = "",
       title = "A. Explicit consent") + 
  scale_y_continuous(expand= c(0,0), limits = c(-0.6, 0.9),
                     breaks = seq(-0.6, 0.9, .3)) +
  coord_flip() +
  theme_minimal() +  
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0, size = 24, face = "italic"),
        plot.subtitle = element_text(hjust = 0, size = 15, face = "italic"),
        plot.title.position = "plot",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(),
        axis.line.x = element_line(size = .3),
        axis.title = element_text(size = 22),
        axis.text = element_text(size = 20),
        axis.text.x = element_text(vjust = -2),
        axis.title.x = element_text(vjust = -3, margin = margin(t = 10, r = 0, b = 0, 
                                                                l = 0)),
        plot.margin = margin(18, 18, 18, 18), 
        aspect.ratio = 0.3) 

# Consent scenario
h2a_y2_lh <- lh_robust(y2_std ~ treatment*age,
                       data = citizens,
                       linear_hypothesis = "treatment + treatment:age = 0")
tidy(h2a_y2_lh)

dt.p <- bind_rows(tidy(h2a_y2_lh) %>% filter(term == "treatment"),
                  tidy(h2a_y2_lh) %>% filter(term == "treatment + treatment:age = 0"))

v.pval <- tidy(h2a_y2_lh) %>% 
  filter(term =="treatment:age") %>% 
  dplyr::select(p.value) %>% 
  pull()

dt.p <- dt.p %>% 
  mutate(diff_p = v.pval,
         xmin = 1,
         xmax = 2, 
         ypos = max(conf.high))

plot_h2a_y2 <- dt.p %>% 
  mutate(label = if_else(term =="treatment", "Old", "Young")) %>% 
  ggplot(aes(y = estimate, x = label, ymin = conf.low,
             ymax = conf.high)) + 
  geom_point(size = 4) +
  geom_errorbar(width = 0.1) +  
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_bracket(data = dt.p %>% slice(1),
               aes(xmin = xmin,
                   xmax = xmax,
                   y.position = ypos,
                   label=if_else(diff_p < .001,
                                 "p < 0.001",
                                 paste0("p = ",(round(diff_p,2))))),
               coord.flip = TRUE,
               label.size = 6,
               color = "#566573",
               tip.length = .02,
               vjust = -.7,
               bracket.nudge.y = 0.05,
               inherit.aes = FALSE) +
  labs(y = "Treatment effect (with 95% CIs)",
       x = "",
       title = "B. Consent scenario") + 
  scale_y_continuous(expand= c(0,0), limits = c(-0.6, 0.9),
                     breaks = seq(-0.6, 0.9, .3)) +
  coord_flip() +
  theme_minimal() +  
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0, size = 24, face = "italic"),
        plot.subtitle = element_text(hjust = 0, size = 15, face = "italic"),
        plot.title.position = "plot",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(),
        axis.line.x = element_line(size = .3),
        axis.title = element_text(size = 22),
        axis.text = element_text(size = 20),
        axis.text.x = element_text(vjust = -2),
        axis.title.x = element_text(vjust = -3, margin = margin(t = 10, r = 0, b = 0, 
                                                               l = 0)),
        plot.margin = margin(18, 18, 18, 18), 
        aspect.ratio = 0.3) 

# Combining plots
grid.arrange(plot_h2a_y1, plot_h2a_y2, nrow = 1)
cate_age <- arrangeGrob(plot_h2a_y1, plot_h2a_y2, nrow = 1, top = "Young people v. old people")
# ggsave(height = 5, width = 20, filename = "./data_analysis/plots/cate_age.png", 
#        cate_age)


# Coefficient plots for gender heterogeneous effects (Table M.3)

# Explicit consent
h2b_y1_lh <- lh_robust(y1_std ~ treatment*gender,
                       data = citizens_young,
                       linear_hypothesis = "treatment + treatment:gender = 0")
tidy(h2b_y1_lh)

dt.p <- bind_rows(tidy(h2b_y1_lh) %>% filter(term=="treatment"),
                  tidy(h2b_y1_lh) %>% filter(term=="treatment + treatment:gender = 0"))

v.pval <- tidy(h2b_y1_lh) %>% 
  filter(term == "treatment:gender") %>% 
  dplyr::select(p.value) %>% 
  pull()

dt.p <- dt.p %>% 
  mutate(diff_p = v.pval,
         xmin = 1,
         xmax = 2, 
         ypos = max(conf.high))

plot_h2b_y1 <- dt.p %>% 
  mutate(label = if_else(term == "treatment", "Man", "Woman")) %>% 
  ggplot(aes(y = estimate, x = label, ymin = conf.low,
             ymax = conf.high)) + 
  geom_point(size = 4) +
  geom_errorbar(width = 0.1) + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_bracket(data=dt.p %>% slice(1),
               aes(xmin = xmin,
                   xmax = xmax,
                   y.position = ypos,
                   label=if_else(diff_p<.001,
                                 "p < 0.001",
                                 paste0("p = ",(round(diff_p,2))))),
               coord.flip = TRUE,
               label.size = 6,
               color="#566573",
               tip.length = .02,
               vjust = -.7,
               bracket.nudge.y = 0.05,
               inherit.aes = FALSE) +
  labs(y = "Treatment effect (with 95% CIs)",
       x = "",
       title = "C. Explicit consent") + 
  scale_y_continuous(expand= c(0,0), limits = c(-0.6, 0.9),
                     breaks = seq(-0.6, 0.9, .3)) +
  coord_flip() +
  theme_minimal() +  
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0, size = 24, face = "italic"),
        plot.subtitle = element_text(hjust = 0, size = 15, face = "italic"),
        plot.title.position = "plot",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(),
        axis.line.x = element_line(size = .3),
        axis.title = element_text(size = 22),
        axis.text = element_text(size = 20),
        axis.text.x = element_text(vjust = -2),
        axis.title.x = element_text(vjust = -3, margin = margin(t = 10, r = 0, b = 0, 
                                                                l = 0)),
        plot.margin = margin(18, 18, 18, 18), 
        aspect.ratio = 0.3) 

# Consent scenario
h2b_y2_lh <- lh_robust(y2_std ~ treatment*gender,
                       data = citizens_young,
                       linear_hypothesis = "treatment + treatment:gender = 0")
tidy(h2b_y2_lh)

dt.p <- bind_rows(tidy(h2b_y2_lh) %>% filter(term == "treatment"),
                  tidy(h2b_y2_lh) %>% filter(term == "treatment + treatment:gender = 0"))

v.pval <- tidy(h2b_y2_lh) %>% 
  filter(term =="treatment:gender") %>% 
  dplyr::select(p.value) %>% 
  pull()

dt.p <- dt.p %>% 
  mutate(diff_p = v.pval,
         xmin = 1,
         xmax = 2, 
         ypos = max(conf.high))

plot_h2b_y2 <- dt.p %>% 
  mutate(label = if_else(term =="treatment", "Man", "Woman")) %>% 
  ggplot(aes(y = estimate, x = label, ymin = conf.low,
             ymax = conf.high)) + 
  geom_point(size = 4) +
  geom_errorbar(width = 0.1) +  
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_bracket(data = dt.p %>% slice(1),
               aes(xmin = xmin,
                   xmax = xmax,
                   y.position = ypos,
                   label=if_else(diff_p < .001,
                                 "p < 0.001",
                                 paste0("p = ",(round(diff_p,2))))),
               coord.flip = TRUE,
               label.size = 6,
               color = "#566573",
               tip.length = .02,
               vjust = -.7,
               bracket.nudge.y = 0.05,
               inherit.aes = FALSE) +
  labs(y = "Treatment effect (with 95% CIs)",
       x = "",
       title = "D. Consent scenario") + 
  scale_y_continuous(expand= c(0,0), limits = c(-0.6, 0.9),
                     breaks = seq(-0.6, 0.9, .3)) +
  coord_flip() +
  theme_minimal() +  
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0, size = 24, face = "italic"),
        plot.subtitle = element_text(hjust = 0, size = 15, face = "italic"),
        plot.title.position = "plot",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(),
        axis.line.x = element_line(size = .3),
        axis.title = element_text(size = 22),
        axis.text = element_text(size = 20),
        axis.text.x = element_text(vjust = -2),
        axis.title.x = element_text(vjust = -3, margin = margin(t = 10, r = 0, b = 0, 
                                                                l = 0)),
        plot.margin = margin(18, 18, 18, 18), 
        aspect.ratio = 0.3) 

# Combining plots
grid.arrange(plot_h2b_y1, plot_h2b_y2, nrow = 1)
cate_gender <- arrangeGrob(plot_h2b_y1, plot_h2b_y2, nrow = 1, top = "Young women vs. young men")
# ggsave(height = 5, width = 20, filename = "./data_analysis/plots/cate_gender.png", 
#        cate_gender)

# Combining coefficient plots on age and gender for Figure 3

# Create a layout matrix with a smaller height for the title row
layout_matrix <- rbind(c(5, 5), c(1, 2), 
                       c(6, 6), c(3, 4))

# Arrange the plots using grid.arrange and the custom layout
p.out <- grid.arrange(plot_h2a_y1, plot_h2a_y2, plot_h2b_y1, plot_h2b_y2, 
                      layout_matrix = layout_matrix,
                      heights = c(1, 4, 1, 4))

# Add a title using grid.text and export
png(filename = "Fig_3_heterogenous_age_gender.png",
    height=600, width=600*1.5)
print(
  grid.arrange(plot_h2a_y1, plot_h2a_y2, plot_h2b_y1, plot_h2b_y2, 
               layout_matrix = layout_matrix,
               heights = c(1, 4, 1, 4))
)

grid.text("Young people vs. old people", x = 0.5, y = 0.925, 
          gp = gpar(fontsize = 24, fontface = "bold"))

grid.text("Young women vs. young men", x = 0.5, y = 0.425, 
          gp = gpar(fontsize = 24, fontface = "bold"))
dev.off()



## Figure 4 and Table M.4: Heterog. eff. for citizens/bureaucrats/pols --------

# citizens <- citizens %>% select(-difference)
administrators_pooled <- administrators %>% select(panel, treatment, y1, y2, 
                                                   age, gender, education, urban, lr) # added code to pool
elected_pooled <- elected %>% select(panel, treatment, y1, y2, 
                                     age, gender, education, urban, lr)  # added code to pool
citizens_pooled <- citizens %>% select(panel, treatment, y1, y2,
                                       age, gender, education, urban, lr)  # added code to pool
pooled <- rbind(administrators_pooled, elected_pooled, citizens_pooled)
pooled$citizens <- ifelse(pooled$panel == "citizens", 1, 0)
pooled$administrators <- ifelse(pooled$panel == "admin", 1, 0)
pooled$elected <- ifelse(pooled$panel == "elected", 1, 0)

pooled <- pooled %>% mutate(panel = factor(panel),
                            panel = relevel(panel, ref = "citizens")) %>%
  dplyr::select(age, gender, education, urban, lr, treatment, 
                panel, administrators, elected, citizens,
                 y1, y2) # added lr, y1, y2

# standardizing outcome variables with a mean of 0 and SD of 1
pooled <- pooled %>% mutate(y1_std = (y1 - mean(y1, na.rm = TRUE))/sd(y1, na.rm = TRUE), 
                                            y2_std = (y2 - mean(y2, na.rm = TRUE))/sd(y2, na.rm = TRUE))

# Checking panel distributions
table(pooled$panel, useNA = "always")


### Table M.4: Heterog. eff. for citizens, politicians, and bureaucrats ###

# Regression models

# Explicit consent
h3_y1 <- lm_robust(y1_std ~ treatment*panel, data = pooled)
summary(h3_y1)

# Consent scenario
h3_y2 <- lm_robust(y2_std ~ treatment*panel, data = pooled)
summary(h3_y2)

# Table output 
screenreg(list(h3_y1, h3_y2), include.ci = FALSE, 
          custom.coef.map = list("(Intercept)" = "Constant", "treatment" = "Treatment", 
                                 "paneladmin" = "Public administrators",
                                 "panelelected" = "Elected representatives",
                                 "treatment:paneladmin" = "Treatment*Public administrators", 
                                 "treatment:panelelected" = "Treatment*Elected representatives"))

texreg(list(h3_y1, h3_y2), 
       include.ci = FALSE, 
       include.rmse = FALSE, 
       file = "Table_M4_heterogenous_citizens_repr_bureau.tex",
       caption = "Exploring heterogeneous treatment effects for public administrators and elected representatives compared with citizens",
       label = "tab:cate_state",
       caption.above = TRUE, 
       custom.header = list("Dependent variable" = 1:2),
       custom.model.name = c("Explicit consent", "Consent scenario"),
       custom.coef.map = list("(Intercept)" = "Constant", "treatment" = "Treatment", 
                              "paneladmin" = "Public administrators",
                              "panelelected" = "Elected representatives",
                              "treatment:paneladmin" = "Treatment*Public administrators", 
                              "treatment:panelelected" = "Treatment*Elected representatives"),
       threeparttable = TRUE,
       custom.note = "\n\\item \\footnotesize \\textit{Note:} OLS regression models with standard errors in parentheses. Sample consists of
       \\item \\footnotesize respondents from the Norwegian Citizen Panel, the Panel of Public Administrators, and  
       \\item \\footnotesize the Panel of Elected Representatives. DVs are standardized with a mean of 0 and a SD
       \\item \\footnotesize of 1. Significant at %stars.")



### Figure 4: Heterog. eff. for citizens, politicians, and bureaucrats ###

# Explicit consent
h3_y1_lh_1 <- lh_robust(y1_std ~ treatment*panel,
                      data = pooled,
                      linear_hypothesis = c("treatment + treatment:paneladmin = 0"))
tidy(h3_y1_lh_1)

h3_y1_lh_2 <- lh_robust(y1_std ~ treatment*panel,
                      data = pooled,
                      linear_hypothesis = c("treatment + treatment:panelelected = 0"))
tidy(h3_y1_lh_2)

dt.p <- bind_rows(tidy(h3_y1_lh_1) %>% filter(term=="treatment"),
                  tidy(h3_y1_lh_1) %>% filter(term=="treatment + treatment:paneladmin = 0"), 
                  tidy(h3_y1_lh_2) %>% filter(term=="treatment + treatment:panelelected = 0"))

v.pval_admin_y1 <- tidy(h3_y1_lh_1) %>% 
  # filter(term == "treatment + treatment:paneladmin = 0") %>% 
  filter(term == "treatment:paneladmin") %>% 
  select(p.value) %>% 
  pull()

v.pval_elected_y1 <- tidy(h3_y1_lh_2) %>% 
  # filter(term == "treatment + treatment:panelelected = 0") %>% 
  filter(term == "treatment:panelelected") %>% 
  select(p.value) %>% 
  pull()

dt.p <- dt.p %>% 
  mutate(# diff_p_admin = v.pval_admin,
         # diff_p_elected = v.pval_elected,
         xmin = 1,
         xmax = 2,
         ypos = max(conf.high))

plot_h3_y1 <- dt.p %>% 
  mutate(label = if_else(term == "treatment", "Citizens", 
                         ifelse(term == "treatment + treatment:paneladmin = 0", "Civil Servants",
                                "Elected Rep.")), 
         label = relevel(factor(label), ref = "Civil Servants")) %>% 
  ggplot(aes(y = estimate, x = label, ymin = conf.low,
             ymax = conf.high)) + 
  geom_point(size = 4) +
  geom_errorbar(width = 0.1) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_bracket(data = dt.p %>% slice(1),
               aes(xmin = xmin,
                   xmax = 1.9,
                   y.position = ypos,
                   label=if_else(v.pval_admin_y1<.001, "p < 0.001", # P-value for interaction term 
                                 if_else(v.pval_admin_y1<.01, "p < 0.01",
                                 paste0("p = ",(round(v.pval_admin_y1,2)))))), # P-value for interaction term
               coord.flip = TRUE,
               label.size = 6,
               color = "#566573",
               tip.length = .02,
               vjust = -0.4,
               bracket.nudge.y = 0.05,
               inherit.aes = FALSE) +
  geom_bracket(data = dt.p %>% slice(3),
               aes(xmin = 2.1,
                   xmax = 3,
                   y.position = ypos,
                   label=if_else(v.pval_elected_y1 <.001, "p < 0.001", 
                                 if_else(v.pval_elected_y1 <.01, "p < 0.01",
                                 paste0("p = ",(round(v.pval_elected_y1, 2)))))),
               coord.flip = TRUE,
               label.size = 6,
               color = "#566573",
               tip.length = .02,
               vjust = -0.4,
               bracket.nudge.y = 0.05,
               inherit.aes = FALSE) +
  labs(y = "Treatment effect (with 95% CIs)",
       x = "",
       title = "A. Explicit consent") + 
  scale_y_continuous(expand= c(0,0), limits = c(-0.3, 0.6),
                     breaks = seq(-0.3, 0.6, .2)) +
  coord_flip() +
  theme_minimal() +  
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0, size = 24, face = "italic"),
        plot.subtitle = element_text(hjust = 0, size = 15, face = "italic"),
        plot.title.position = "plot",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(),
        axis.line.x = element_line(linewidth = .3),
        axis.title = element_text(size = 22),
        axis.text = element_text(size = 20),
        axis.text.x = element_text(vjust = -2),
        axis.title.x = element_text(vjust = -3, margin = margin(t = 10, r = 0, b = 0, 
                                                                l = 0)),
        plot.margin = margin(18, 18, 18, 18))

# Consent scenario
h3_y2_lh_1 <- lh_robust(y2_std ~ treatment*panel,
                      data = pooled,
                      linear_hypothesis = c("treatment + treatment:paneladmin = 0"))
tidy(h3_y2_lh_1)

h3_y2_lh_2 <- lh_robust(y2_std ~ treatment*panel,
                      data = pooled,
                      linear_hypothesis = c("treatment + treatment:panelelected = 0"))
tidy(h3_y2_lh_2)

dt.p <- bind_rows(tidy(h3_y2_lh_1) %>% filter(term=="treatment"),
                  tidy(h3_y2_lh_1) %>% filter(term=="treatment + treatment:paneladmin = 0"), 
                  tidy(h3_y2_lh_2) %>% filter(term=="treatment + treatment:panelelected = 0"))

v.pval_admin_y2 <- tidy(h3_y2_lh_1) %>% 
  # filter(term == "treatment + treatment:paneladmin = 0") %>% 
  filter(term == "treatment:paneladmin") %>% 
  select(p.value) %>% 
  pull()

v.pval_elected_y2 <- tidy(h3_y2_lh_2) %>% 
  # filter(term == "treatment + treatment:panelelected = 0") %>% 
  filter(term == "treatment:panelelected") %>% 
  select(p.value) %>% 
  pull()

dt.p <- dt.p %>% 
  mutate(# diff_p_admin = v.pval_admin,
         # diff_p_elected = v.pval_elected,
         xmin = 1,
         xmax = 2,
         ypos = max(conf.high))

plot_h3_y2 <- dt.p %>% 
  mutate(label = if_else(term == "treatment", "Citizens", 
                         ifelse(term == "treatment + treatment:paneladmin = 0", "Civil Servants",
                                "Elected Rep.")), 
         label = relevel(factor(label), ref = "Civil Servants")) %>% 
  ggplot(aes(y = estimate, x = label, ymin = conf.low,
             ymax = conf.high)) + 
  geom_point(size = 4) +
  geom_errorbar(width = 0.1) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_bracket(data = dt.p %>% slice(1),
               aes(xmin = xmin,
                   xmax = 1.9,
                   y.position = ypos,
                   label=if_else(v.pval_admin_y2<.001, "p < 0.001",
                                 if_else(v.pval_admin_y2<.01, "p < 0.01",
                                 paste0("p = ",(round(v.pval_admin_y2,2)))))),
               coord.flip = TRUE,
               label.size = 6,
               color = "#566573",
               tip.length = .02,
               vjust = -0.4,
               bracket.nudge.y = 0.05,
               inherit.aes = FALSE) +
  geom_bracket(data = dt.p %>% slice(3),
               aes(xmin = 2.1,
                   xmax = 3,
                   y.position = ypos,
                   label=if_else(v.pval_elected_y2<.001, "p < 0.001",
                                 if_else(v.pval_elected_y2<.01, "p < 0.01",
                                 paste0("p = ",(round(v.pval_elected_y2,2)))))),
               coord.flip = TRUE,
               label.size = 6,
               color = "#566573",
               tip.length = .02,
               vjust = -0.4,
               bracket.nudge.y = 0.05,
               inherit.aes = FALSE) +
  labs(y = "Treatment effect (with 95% CIs)",
       x = "",
       title = "B. Consent scenario") + 
  scale_y_continuous(expand= c(0,0), limits = c(-0.3, 0.6),
                     breaks = seq(-0.3, 0.6, .2)) +
  coord_flip() +
  theme_minimal() +  
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0, size = 24, face = "italic"),
        plot.subtitle = element_text(hjust = 0, size = 15, face = "italic"),
        plot.title.position = "plot",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(),
        axis.line.x = element_line(linewidth = .3),
        axis.title = element_text(size = 22),
        axis.text = element_text(size = 20),
        axis.text.x = element_text(vjust = -2),
        axis.title.x = element_text(vjust = -3, margin = margin(t = 10, r = 0, b = 0, 
                                                                l = 0)),
        plot.margin = margin(18, 18, 18, 18))

# Combining plots
grid.arrange(plot_h3_y1, plot_h3_y2, nrow = 1)
cate_state <- arrangeGrob(plot_h3_y1, plot_h3_y2, nrow = 1)
ggsave(height = 5, width = 15, filename = "Fig_4_heterogenous_citizens_repr_bureau.png", 
       cate_state)



### Further analysis ----

# Proportions for baseline support
prop.table(table(citizens_baseline$y1))
prop.table(table(citizens_baseline$y2))

# Means for control v. treatment group
mean(citizens$y1[citizens$treatment == 0], na.rm = TRUE)
mean(citizens$y1[citizens$treatment == 1], na.rm = TRUE)

mean(citizens$y2[citizens$treatment == 0], na.rm = TRUE)
mean(citizens$y2[citizens$treatment == 1], na.rm = TRUE)

#### Power analysis #### 
set.seed(9876)

# Testing power: Citizen Panel
pwr.t.test(n = 750, sig.level = 0.05, power = 0.90, type = "two.sample") # 0.17
pwr.t.test(n = 750, sig.level = 0.025, power = 0.90, type = "two.sample") # 0.18

d.05 <- pwr.t.test(n = 750, sig.level = 0.05, power = 0.90, type = "two.sample")$d
d.025 <- pwr.t.test(n = 750, sig.level = 0.025, power = 0.90, type = "two.sample")$d

# Setting power 
power <- 0.90

# Creating function for power plot
power.curve.citizen <- function(n, effect){
  samp.out <- NULL
  
  for(i in 1:length(n)){
    
    power <-  pwr.t.test(d=effect[i], 
                         n=n[i]/2, #NB
                         sig.level=.05,
                         type="two.sample")$power
    power <-  data.frame(effect.size=effect[i],sample.size=n[i],power=power)
    samp.out <- rbind(samp.out, power)
  }
  
  return(samp.out)
}

# Setting sample size and effect sizes
effect.05 <- rep(sort(c(d.05, seq(0, 0.4, 0.1))), 3) #Vector of MDEs

samples <- lapply(list(1300, 1500, 1700), rep, times = 6, simplify = TRUE)
samples <- unlist(samples)

# Plotting power analysis
citizen_power <- power.curve.citizen(samples, effect.05)
citizen_power$sample.size <- as.factor(citizen_power$sample.size)
citizen_power_plot <- ggplot(citizen_power,
                             aes(effect.size, power, 
                                 linetype = sample.size,
                                 shape = sample.size,
                                 color = sample.size)) + 
  geom_hline(yintercept = 0.9, linetype = "dotted", color="darkgrey") + 
  geom_vline(xintercept = d.05, linetype = "dotted", color="darkgrey") +
  annotate("text", x = d.05+.005, y = 0, label = paste0("MDE =\n",round(d.05,2)),
           vjust = .2, hjust = 0, size = 4, color="darkgrey") +
  annotate("text", x = 0, y = 0.9, label = paste0("Power = 0.9"),
           vjust = -.5, hjust = .15, size = 4, color="darkgrey") +
  geom_line(linewidth = .5) + 
  geom_point(size = 2) + 
  scale_linetype_manual(values=c("twodash","solid","longdash")) + 
  scale_color_manual(values=c("darkgrey","black","darkgrey")) + 
  scale_shape_manual(values=c(17,16,15)) + 
  scale_y_continuous(breaks=c(seq(0,1,.25)), limits=c(0,1)) + 
  theme(axis.title.x = element_text(size = 12),
        axis.title.y = element_text(size = 12),
        axis.text = element_text(size = 12),
        plot.title = element_text(hjust=0, size=16, face="italic"),
        plot.title.position = "plot",
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        legend.position = "none", 
        plot.margin = margin(0, 0, 0, 0, "pt"))  +
  labs(x = "Minimal detectable effect (MDE)", y = "Power", 
       linetype = "Sample size", 
       title = "A. Alpha = 0.05")

ggsave(height = 5, width = 6, filename = "citizen_power_plot.png", 
       citizen_power_plot)

# Creating function for power plot (penalized)
power.curve.citizen.penalized <- function(n, effect){
  samp.out <- NULL
  
  for(i in 1:length(n)){
    
    power <-  pwr.t.test(d=effect[i], 
                         n=n[i]/2, #NB
                         sig.level=.025,
                         type="two.sample")$power
    power <-  data.frame(effect.size=effect[i],sample.size=n[i],power=power)
    samp.out <- rbind(samp.out, power)
  }
  
  return(samp.out)
}

# Setting sample size 
effect.025 <- rep(sort(c(d.025, seq(0, 0.4, 0.1))), 3) #Vector of MDEs

samples <- lapply(list(1300, 1500, 1700), rep, times = 6, simplify = TRUE)
samples <- unlist(samples)

# Plotting power analysis
citizen_power_penalized <- power.curve.citizen.penalized(samples, effect.025)
citizen_power_penalized$sample.size <- as.factor(citizen_power_penalized$sample.size)
citizen_power_plot_penalized <-
  ggplot(citizen_power_penalized,
         aes(effect.size, power, 
             linetype = sample.size,
             shape = sample.size,
             color = sample.size)) + 
  geom_hline(yintercept = 0.9, linetype = "dotted", color="darkgrey") + 
  geom_vline(xintercept = d.025, linetype = "dotted", color="darkgrey") +
  annotate("text", x = d.025+.005, y = 0, label = paste0("MDE =\n",round(d.025,2)),
           vjust = .2, hjust = 0, size = 4, color="darkgrey") +
  annotate("text", x = 0, y = 0.9, label = paste0("Power = 0.9"),
           vjust = -.5, hjust = .15, size = 4, color="darkgrey") +
  geom_line(linewidth = .5) + 
  geom_point(size = 2) + 
  scale_linetype_manual(values=c("twodash","solid","longdash"),
                        labels=c("1,300", "1,500 (exp.)", "1,700"),
                        name="Sample size") + 
  scale_color_manual(values=c("darkgrey","black","darkgrey"),
                     labels=c("1,300", "1,500 (exp.)", "1,700"),
                     name="Sample size") + 
  scale_shape_manual(values=c(17,16,15),
                     labels=c("1,300", "1,500 (exp.)", "1,700"),
                     name="Sample size") +   
  scale_y_continuous(breaks=c(seq(0,1,.25)), limits=c(0,1)) + 
  guides(fill = guide_legend(ncol=1),
         color = guide_legend(ncol=1),
         shape = guide_legend(ncol=1),
         linetype = guide_legend(ncol=1)) +
  theme(axis.title.x = element_text(size = 12),
        axis.title.y = element_text(size = 12),
        axis.text = element_text(size = 12),
        plot.title = element_text(hjust=0, size=16, face="italic"),
        plot.title.position = "plot",
        legend.text = element_text(size = 9),
        # legend.title = element_text(size = 9),
        # legend.key.spacing.y = unit(.1, 'pt'),
        legend.box = "vertical",
        legend.position = "inside",
        legend.position.inside = c(0.83, 0.5),
        legend.background = element_rect(color="black", linewidth=.1),
        legend.title = element_text(margin=margin(b=.1), size=9),
        legend.key.spacing.y = unit(-7, "pt"),
        plot.margin = margin(0, 0, 0, 0, "pt"))  +
  # scale_linetype_manual(values = c("1300" = 4, "1500" = 1, "1700" = 3), 
  #                       labels = c("1300", "1500 (Exp.)", "1700")) +
  labs(x = "Minimal detectable effect (MDE)", y = "Power", 
       # linetype = "Sample size", 
       title = "B. Alpha = 0.025")

ggsave(height = 5, width = 6, filename = "citizen_power_plot_penalized.png", 
       citizen_power_plot_penalized)

# Plot for paper

grid.arrange(citizen_power_plot, citizen_power_plot_penalized, nrow = 1)
citizen_power_bh <- arrangeGrob(citizen_power_plot, citizen_power_plot_penalized, nrow = 1)
ggsave(height = 3, width = 8, filename = "citizen_power_bh.png", 
       citizen_power_bh)